An Reproducible Example of Colocalization Analysis
Yonghui Dong
(1) Load packages and read data
The dataset is downloaded from Metabolights. Accession number: MTBLS487. The data was obtained in positive ion mode over an m/z range of 300–2000 using a Thermo MALDI LTQ Orbitrap XL instrument (Thermo Fisher Scientific, Bremen, Germany). The spatial resolution was 50 micrometers. For more information about the dataset, readers are kindly referred to the original publication (Hall et al., 2017).
## Main package and the version, R version 4.2.0
library(Cardinal) #(v2.14.0)
library(BiocParallel) #(v1.30.4)
library(plotly) #(v4.10.1)
## Read imzML dataset
path <- paste0("../Data/test_POS", ".imzML")
Brain <- readMSIData(path, resolution = 3, units = "ppm",
mass.range = c(300, 2000), attach.only = T)(2) Pre-processing
Brain2 <-
Brain %>%
normalize(method = "tic") %>%
peakPick(method = "mad", SNR = 6) %>%
peakAlign(tolerance = 3, units = "ppm") %>%
process(BPPARAM = SerialParam())Check the ion image of m/z 844.527 before and after pre-processing in order to make sure that chosen pre-processing parameters are fine.
darkmode()
##(2.1) before pre-processing
image(Brain, mz = 844.527, smooth.image = "gaussian", plusminus = 0.03, colorscale = magma,
contrast.enhance = "suppression", normalize.image = "linear")
##(2.2) after pre-processing
image(Brain2, mz = 844.527, smooth.image = "gaussian", plusminus = 0.03, colorscale = magma,
contrast.enhance = "suppression", normalize.image = "linear")
(3) Co-localization
Here first 100 colocalized ions are retained. Three colocalization
measures were calculated. correlation is Pearson’s
correlation, M1 is match score, and M2 is
Manders’ colocalization coefficients. Here Pearson’s correlation is
used.
#(3.1) co-localization PC36:1
coloc_826 <- colocalized(Brain2, mz= 826.572, n = 100, BPPARAM = SerialParam())
#(3.2) co-localization PC38:6
coloc_844 <- colocalized(Brain2, mz= 844.525, n = 100, BPPARAM = SerialParam())
#(3.3) co-localization PC40:6
coloc_872 <- colocalized(Brain2, mz= 872.557, n = 100, BPPARAM = SerialParam())(4) Extract in-source fragments, and create pseudo MS/MS spectrum
#(4.1) PC36:1
mycoloc <- as.data.frame(coloc_826)
MSe <- mycoloc[mycoloc$correlation >= 0.90,]
## extract intensity
int = as.vector(rep(NA, dim(MSe)[1]))
for (i in 1:length(MSe$mz)) {
int[i] <- sum(spectra(Brain2)[features(Brain2, mz = MSe$mz[i]),])
}
## interactive plot
spec = cbind.data.frame(mz = MSe$mz, Int = int)
p = ggplot(spec, aes(x = mz, ymax = Int/max(Int)*100, ymin = 0, colour = "red")) +
geom_linerange() +
scale_y_continuous(expand = c(0, 0), limits = c(0, 100*1.1)) +
theme_bw() +
ggtitle("m/z 826.577")
ggplotly(p)#(4.2) PC38:6
mycoloc <- as.data.frame(coloc_844)
MSe <- mycoloc[mycoloc$correlation >= 0.90,]
## extract intensity
int = as.vector(rep(NA, dim(MSe)[1]))
for (i in 1:length(MSe$mz)) {
int[i] <- sum(spectra(Brain2)[features(Brain2, mz = MSe$mz[i]),])
}
## interactive plot
spec = cbind.data.frame(mz = MSe$mz, Int = int)
p = ggplot(spec, aes(x = mz, ymax = Int/max(Int)*100, ymin = 0, colour = "red")) +
geom_linerange() +
scale_y_continuous(expand = c(0, 0), limits = c(0, 100*1.1)) +
theme_bw() +
ggtitle("m/z 844.528")
ggplotly(p)#(4.3) PC40:6
mycoloc <- as.data.frame(coloc_872)
MSe <- mycoloc[mycoloc$correlation >= 0.90,]
## extract intensity
int = as.vector(rep(NA, dim(MSe)[1]))
for (i in 1:length(MSe$mz)) {
int[i] <- sum(spectra(Brain2)[features(Brain2, mz = MSe$mz[i]),])
}
## interactive plot
spec = cbind.data.frame(mz = MSe$mz, Int = int)
p = ggplot(spec, aes(x = mz, ymax = Int/max(Int)*100, ymin = 0, colour = "red")) +
geom_linerange() +
scale_y_continuous(expand = c(0, 0), limits = c(0, 100*1.1)) +
theme_bw()+
ggtitle("m/z 872.559")
ggplotly(p)(5) Plot MALDI images
#(5.1) PC36:1
mycoloc <- as.data.frame(coloc_826)
MSe <- mycoloc[mycoloc$correlation >= 0.90,]
pdf(file = file.path("Result/colocalization_826.pdf"), onefile = TRUE)
for(i in 1:dim(MSe)[1]){
darkmode()
print(image(Brain2, mz = MSe$mz[i], smooth.image = "gaussian",
plusminus = 0.003, colorscale=magma,
contrast.enhance="suppression", normalize.image = "linear"))
legend("topleft", legend= paste("correlation = ", round(MSe$correlation[i], 2)))
}
dev.off()
#(5.2) PC38:6
mycoloc <- as.data.frame(coloc_844)
MSe <- mycoloc[mycoloc$correlation >= 0.90,]
pdf(file = file.path("Result/colocalization_844.pdf"), onefile = TRUE)
for(i in 1:dim(MSe)[1]){
darkmode()
print(image(Brain2, mz = MSe$mz[i], smooth.image = "gaussian",
plusminus = 0.003, colorscale=magma,
contrast.enhance="suppression", normalize.image = "linear"))
legend("topleft", legend= paste("correlation = ", round(MSe$correlation[i], 2)))
}
dev.off()
#(5.3) PC38:6
mycoloc <- as.data.frame(coloc_872)
MSe <- mycoloc[mycoloc$correlation >= 0.90,]
pdf(file = file.path("Result/colocalization_872.pdf"), onefile = TRUE)
for(i in 1:dim(MSe)[1]){
darkmode()
print(image(Brain2, mz = MSe$mz[i], smooth.image = "gaussian",
plusminus = 0.003, colorscale=magma,
contrast.enhance="suppression", normalize.image = "linear"))
legend("topleft", legend= paste("correlation = ", round(MSe$correlation[i], 2)))
}
dev.off()Session information
sessioninfo::session_info()## ─ Session info ───────────────────────────────────────────────────────────────
## setting value
## version R version 4.2.0 (2022-04-22)
## os macOS Big Sur/Monterey 10.16
## system x86_64, darwin17.0
## ui X11
## language (EN)
## collate en_US.UTF-8
## ctype en_US.UTF-8
## tz Asia/Jerusalem
## date 2022-12-10
## pandoc 2.19.2 @ /Applications/RStudio.app/Contents/MacOS/quarto/bin/tools/ (via rmarkdown)
##
## ─ Packages ───────────────────────────────────────────────────────────────────
## package * version date (UTC) lib source
## abind 1.4-5 2016-07-21 [1] CRAN (R 4.2.0)
## assertthat 0.2.1 2019-03-21 [1] CRAN (R 4.2.0)
## biglm 0.9-2.1 2020-11-27 [1] CRAN (R 4.2.0)
## Biobase 2.56.0 2022-04-26 [1] Bioconductor
## BiocGenerics * 0.42.0 2022-04-26 [1] Bioconductor
## BiocParallel * 1.30.4 2022-10-13 [1] Bioconductor
## bitops 1.0-7 2021-04-24 [1] CRAN (R 4.2.0)
## bslib 0.4.1 2022-11-02 [1] CRAN (R 4.2.0)
## cachem 1.0.6 2021-08-19 [1] CRAN (R 4.2.0)
## Cardinal * 2.14.0 2022-04-26 [1] Bioconductor
## cli 3.4.1 2022-09-23 [1] CRAN (R 4.2.0)
## codetools 0.2-18 2020-11-04 [1] CRAN (R 4.2.0)
## colorspace 2.0-3 2022-02-21 [1] CRAN (R 4.2.0)
## crosstalk 1.2.0 2021-11-04 [1] CRAN (R 4.2.0)
## data.table 1.14.6 2022-11-16 [1] CRAN (R 4.2.0)
## DBI 1.1.3 2022-06-18 [1] CRAN (R 4.2.0)
## digest 0.6.30 2022-10-18 [1] CRAN (R 4.2.0)
## dplyr 1.0.10 2022-09-01 [1] CRAN (R 4.2.0)
## EBImage * 4.38.0 2022-05-15 [1] Bioconductor
## ellipsis 0.3.2 2021-04-29 [1] CRAN (R 4.2.0)
## evaluate 0.18 2022-11-07 [1] CRAN (R 4.2.0)
## fansi 1.0.3 2022-03-24 [1] CRAN (R 4.2.0)
## farver 2.1.1 2022-07-06 [1] CRAN (R 4.2.0)
## fastmap 1.1.0 2021-01-25 [1] CRAN (R 4.2.0)
## fftwtools 0.9-11 2021-03-01 [1] CRAN (R 4.2.0)
## generics 0.1.3 2022-07-05 [1] CRAN (R 4.2.0)
## ggplot2 * 3.3.6 2022-05-03 [1] CRAN (R 4.2.0)
## glue 1.6.2 2022-02-24 [1] CRAN (R 4.2.0)
## gtable 0.3.1 2022-09-01 [1] CRAN (R 4.2.0)
## highr 0.9 2021-04-16 [1] CRAN (R 4.2.0)
## htmltools 0.5.3 2022-07-18 [1] CRAN (R 4.2.0)
## htmlwidgets 1.5.4 2021-09-08 [1] CRAN (R 4.2.0)
## httr 1.4.4 2022-08-17 [1] CRAN (R 4.2.0)
## irlba 2.3.5.1 2022-10-03 [1] CRAN (R 4.2.0)
## jpeg 0.1-9 2021-07-24 [1] CRAN (R 4.2.0)
## jquerylib 0.1.4 2021-04-26 [1] CRAN (R 4.2.0)
## jsonlite 1.8.3 2022-10-21 [1] CRAN (R 4.2.0)
## knitr 1.41 2022-11-18 [1] CRAN (R 4.2.0)
## labeling 0.4.2 2020-10-20 [1] CRAN (R 4.2.0)
## lattice 0.20-45 2021-09-22 [1] CRAN (R 4.2.0)
## lazyeval 0.2.2 2019-03-15 [1] CRAN (R 4.2.0)
## lifecycle 1.0.3 2022-10-07 [1] CRAN (R 4.2.0)
## locfit 1.5-9.6 2022-07-11 [1] CRAN (R 4.2.0)
## magrittr 2.0.3 2022-03-30 [1] CRAN (R 4.2.0)
## MASS 7.3-58.1 2022-08-03 [1] CRAN (R 4.2.0)
## Matrix 1.5-3 2022-11-11 [1] CRAN (R 4.2.0)
## matter 1.22.0 2022-04-26 [1] Bioconductor
## mclust 6.0.0 2022-10-31 [1] CRAN (R 4.2.0)
## munsell 0.5.0 2018-06-12 [1] CRAN (R 4.2.0)
## nlme 3.1-160 2022-10-10 [1] CRAN (R 4.2.0)
## pillar 1.8.1 2022-08-19 [1] CRAN (R 4.2.0)
## pkgconfig 2.0.3 2019-09-22 [1] CRAN (R 4.2.0)
## plotly * 4.10.1 2022-11-07 [1] CRAN (R 4.2.0)
## png 0.1-7 2013-12-03 [1] CRAN (R 4.2.0)
## ProtGenerics * 1.28.0 2022-04-26 [1] Bioconductor
## purrr 0.3.5 2022-10-06 [1] CRAN (R 4.2.0)
## R6 2.5.1 2021-08-19 [1] CRAN (R 4.2.0)
## RCurl 1.98-1.9 2022-10-03 [1] CRAN (R 4.2.0)
## rlang 1.0.6 2022-09-24 [1] CRAN (R 4.2.0)
## rmarkdown 2.18 2022-11-09 [1] CRAN (R 4.2.0)
## rstudioapi 0.14 2022-08-22 [1] CRAN (R 4.2.0)
## S4Vectors * 0.34.0 2022-04-26 [1] Bioconductor
## sass 0.4.4 2022-11-24 [1] CRAN (R 4.2.0)
## scales 1.2.1 2022-08-20 [1] CRAN (R 4.2.0)
## sessioninfo 1.2.2 2021-12-06 [1] CRAN (R 4.2.0)
## signal 0.7-7 2021-05-25 [1] CRAN (R 4.2.0)
## sp 1.5-1 2022-11-07 [1] CRAN (R 4.2.0)
## stringi 1.7.8 2022-07-11 [1] CRAN (R 4.2.0)
## stringr 1.4.1 2022-08-20 [1] CRAN (R 4.2.0)
## tibble 3.1.8 2022-07-22 [1] CRAN (R 4.2.0)
## tidyr 1.2.1 2022-09-08 [1] CRAN (R 4.2.0)
## tidyselect 1.2.0 2022-10-10 [1] CRAN (R 4.2.0)
## tiff 0.1-11 2022-01-31 [1] CRAN (R 4.2.0)
## utf8 1.2.2 2021-07-24 [1] CRAN (R 4.2.0)
## vctrs 0.5.1 2022-11-16 [1] CRAN (R 4.2.0)
## viridisLite 0.4.1 2022-08-22 [1] CRAN (R 4.2.0)
## withr 2.5.0 2022-03-03 [1] CRAN (R 4.2.0)
## xfun 0.35 2022-11-16 [1] CRAN (R 4.2.0)
## yaml 2.3.6 2022-10-18 [1] CRAN (R 4.2.0)
##
## [1] /Library/Frameworks/R.framework/Versions/4.2/Resources/library
##
## ──────────────────────────────────────────────────────────────────────────────